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ABSTRACT 


This  paper  presents  an  accuracy  analysis  of  a  suggested  approximate 
confidence  interval  for  system  maintainability  parameters. 

Technically,  the  simulation  demonstrates  feasible  ranges  of"  parameter 
applicability  for  a  fit  of  linear  combinations  of  generated  gamma  variates 
to  the  gamma  distribution,  using  the  method  of  moments. 

The  simulation  has  application  to  the  classical  confidence-  interval 
for  mean  time  to  repair  of  a  series  system,  under  the  assumptions  of 
gamma  distributed  repair  times,  and  method  of  moments  estimators. 

The  paper  provides  no  validated  conclusions  although  it  does  display 
parameters  and  ranges  of  apparent  extremely  high  model  validity.. 
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I.   INTRODUCTION 

The  purpose  of  this  paper  is  to  present  a  computer  simulation,  in 
order  to  investigate  the  adequacy  of  a  suggested  approximate  confidence 
interval  for  system  maintainability  parameters.   In  a  technical  sense, 
it  will  demonstrate  feasible  ranges  of  parameter  applicability  for  a  fit 
of  linear  combinations  of  generated  gamma  variates  to  the  gamma  distribu- 
tion, using  the  method  of  moments.   The  application  of  the  simulation  is 
to  the  classical  confidence  interval  for  mean  time  to  repair  a  system, 
under  the  assumptions  of  series  components,  gamma  distributed  repair 
times,  and  method  of  moments  estimators.   The  parameters  of  the  distribu- 
tions that  were  investigated  were  selected  so  that  the  maintainability 
curves  would  be  representative  of  those  curves  presently  of  interest  in 
industry. 


II.   SUMMARY 

A  procedure  for  a  fit  of  linear  combinations  of  generated  gammas 
variates  to  the  gamma  distribution,  using  the  method  of  moments  is: 
presented  in  Chapter  III.   A  simulation  of  this  procedure  is  proposed 
in  Chapter  IV  and  the  results  are  tabulated  in  Tables  I  and  II  of 
Chapter  V. 


III.   THE  STATISTICAL  MODEL 

A.  SERIES  ASSUMPTION 

Suppose  we  are  given  a  system  of  components  which  are  not  necessarily 
series,  where  the  time  to  failure  of  the  ith  component  of  the  system  has  the 
exponential  (X.)  distribution. 

Then, 

1)  for  the  case  where  there  are  k.  identical  components  in.  series  of 
type  i  within  the  system,  set  X.'  =  k.X.. 

2)  for  the  case  where  there  are  h.  identical  components  in  parallel 

J 

of  type  j  within  the  system,  set  X.'  =  X.  j 

Therefore,  for  some  parallel-series  combination  there  is  obtainable 
a  series  system  of  uniquely  different  components  whose  failure  ratios  are 
assumed  to  be  exponential. 

Thus,  the  system  being  described  is  a  series  system  with  k  different 

components  and  failure  rates  X.,  X„,  •••,  X.fT  X  %  ----.  X,  where 

12xj         k 

all  the  X.  are  in  the  same  units. 

1 

B.  FAILURE  RATE  ASSUMPTION 

Let  X.  denote  the  time  to  failure  of  component  i,  where  X.  has  the 
exponential  (X.)  distribution.   Thus, 

fv  (x. ;  X.)  =  X.e"Xixi  (3-1) 

X.    11      l 
l 

We  further  assume  that  the  system  fails  when  exactly  one  component 

fails  or  the  component  which  caused  the  system  failure  takes  the  longest 

time  to  repair. 


C.   GAMMA  ASSUMPTION 

Let  T.  denote  the  time  required  to  repair  the  ith  component  and 
suppose  that  T.  has  the  gamma  distribution.   Thus,  T.  **  T    (t.  ;  a.  ,.  tr.).. 


Therefore, 


b.   b.-l  -a.t. 
a.  1  t .  1   e  11 

f   (t.;  a.,  b.)  -  -i i^ (i-2) 


is  the  density  equation  with 


b. 
E[T  ]  =  -1  =  y,  C3H3) 

X       ci  ,        1 

i 
which  is  the  mean  time  to  repair  (MTTR)  for  the  ith  component. 

D.  SYSTEM  MTTR  OBTAINED 

Let  9  denote  the  mean  time  to  repair  for  the  system.   Therefore, 

9  =  E  P  [component  i  fails  first]  y.  ('3-4) 

i  X 

and,  if  the  assumption  of  a  series  system  is  valid,  then 

k 
X  =   E   X. 
1-1   * 

and 

k  X.  b. 

9  =  MTTR  =   E  (~)  (— )  (3-5) 

i-1  X  ai 

E.  CONFIDENCE  LIMIT  OBTAINED 

An  estimated  upper  confidence  limit  for  system  maintainability,  denoted 

by  9  ,  has  been  proposed.   9  was  derived  by  the  following  procedure: 

For  any  component  i,  T,_ ,  T_ _,  T, _,  •••  T.    is  a  random  sample  on  T., 

11    12    U       m.  i 

i 

the  time  to  repair  component  i.   Failure  data  is  also  available  on  component 

i,  so  that  an  estimate  X.  of  X.  is  possible. 

li 


Thus,  9  can  be  estimated  by  9  where 


k  X.  _ 

&  -   E   71  T  (3r-6) 

i=l  X         X 


and 


T,  -  -IT.. 
1    n.  .   it 
i  J 

Moreover,  because  of  the  large  amount  of  industrial  testing,,  it ~ can 

be  assumed  that  X.  =  X.,  and  for  purposes  of  derivation  we  shall,  treat 
11 

the  X.  as  though  they  are  constants.  Thus, 

k   X. 
9  =  £  T^T  (3-7) 

i=l  A   1 


k  X.  b. 

E  [9]  =  Z  -±   (— )  (3-8) 

.  -»  X   a. 
i=l     i 

and 


k  X.  2  b. 
var  (9)  =  I  (^)  (-^-)  .  (3-9) 

i=l      a. 

i 

We  shall  fit  a  two  parameter  gamma  to  the  density  of  9  and  obtain 
the  upper  confidence  limit  from  this  fitted  distribution.   Thus  we  are 
assuming  9   T(a,b)  and  use  will  be  made  of  the  method  of  moments  to  make 
the  fit.   Through  use  of  formulas  (3-8)  and  (3-9) 


k  X.  b. 
E  [9]  =  Z   T-i  (-%  =  7  (3-10) 

1=1    i 


and 


k   X.  2  b.. 
var  (9)  =  Z      (-^      (-^>  =  ^y  (0341) 

i=l        a .     a 

l 
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The  parameters  a  and  b  can  be  solved  for,  as  follows 


E   (~)  (j1)  <MZ) 


a  = 


k   X.   b 

*  (3T>  <I. 

i=l i_ 

k   X.2  b. 

i=l       a. 

1 


and 


k  X.  b.   2 

b  =  ^  X.lb.  ■     «M3J 

i=l       a. 

i 

For  the  available  data,  this  then  becomes: 


k  X. 
E  v1  T. 

b  =  1=1  \   = (3-44) 

k  X. 

E  fei)   S.2 

i=l  X     x 

where , 

n. 

1    X       -  2 
S   =  -t-r-   E   (T,  -  T)Z 

1    ni_1  j=l    J 

2 
If  b  is  an  integer,  29a  has  the  y,_,  n  distribution:  therefore,,  define 

C2b) 

(2b)  as  the  integer  closest  to  2b  and  the  following  approximate  confidence 
interval  can  be  formed: 


10 


l-a  =  P  [20a  >  X2(1^)(2b)l  C*-15) 


2 
=  P  r29a  >  xd-ct)(2b) 

l(2b)      (2b)    J 


.  P  C1L  <  ^1 ] 

a    2 

X(l-a)  (2b) 


*  P  E|  <  92(2b)   :  ]  CK6) 

X(l-a)  (2b) 

where  the  approximate  equality  compensates  for  the  use  of  estimators ► 

Finally,  the  choice  of  a  =  .20  results  in  the  following: 


.80  *P  [±<   -f-^ ]•  CW7) 

a    i 


X(,80)(2b) 


Thus,  the  desired  suggested  approximate  80%  confidence  limit  for  system 

2 
maintainability  is  9(2b)/x,  0wov\  • 

{ .o)  (,/b; 
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IV,   THE  SIMULATION  PROCEDURE 

As  explained  in  Chapter  III,  the  system  to  be  simulated  consists  of 
k  components  in  logical  series  with  system  maintainability  (9)  expressed 
as : 


k   X.   b. 
9  =  Z      (ri)  (-±) 


i=l 


a. 

1 


(^-1) 


where  b./a.  is  the  true  maintainability  of  the  ith  component  of  the 
xi  . 

system. 

Denote  an  upper  confidence  limit  for  0  by  9  .   If  &  is  in  fact  the 


u 


exact  100  (l-a)%  upper  confidence  limit  for  9,  then 


b   9  (2b) 
La    2 


•]  =  1  -  a 


A(l-a)(2b) 

holds  to  a  reasonably  close  approximation. 

In  fact,  b/a  should  then  be  the  ath  percentile  point  of  the  simulated 

distribution  of  9  . 
u 


b 
a 


The  choice  of  a  =  .20  defined  9  /oN  as  the  20th  percentile  point  of 

u(2) 


the  distribution  of  9  . 

u 
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In  order  to  investigate  this  distribution,  a  computer  was  utilized 

to  generate  the  required  gamma  variates.   Then,  500  values  of  9  were 

u 

computer  and  ordered  such  that 

Kl  *  Kl  *   »u3  <  •  • '  «  »u500  '      '**> 

Since  it  was  desired  to  display  the  80%  upper  confidence  limit  for 

8  (0  /9\)>  which  implies  80%  of  the  ordered  values  be  greater  than  0  ,..,, 

the  20th  percentile  point  of  the  9  distribution  was  found.  This  ZOth 

u 

percentile  point  is  the  100th  ordered  value  in  formula  (4-2)  above  and. 

if  the  procedure  is  correct,  should  approximately  equal  — . 

a 

The  primary  measure  of  accuracy  for  the  simulation  will  be  the  value 
of 

Ul0°   a'  ((4-3) 


b/a 

which  is  an  expression  relating  the  estimated  value  of  system  maintain- 
ability  (9   _  )  to  the  gamma  value  for  MTTR  (b/a).   Thus,  the  accuracy 
of  the  simulation  is  presented  in  the  notion  of  relative  error. 
Also,  the  statistic 


[//9  .  >  b/a] 

H2_= (4-4) 

500 

will  be  computed  in  order  to  display  the  relative  error  between  the  true 
value  of  system  maintainability  (b/a)  and  the  number  of  generated  esti- 
mates of  this  value  (9  .). 

The  analysis  for  this  proposed  method  was  conducted  using  different 
combinations  of  values  of  the  gamma  input  parameters,  varying  A.  values, 
one  of  two  values  for  k,  and  n.  values  of  10,  20,  50,  and  100. 
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A.   THE  SIMULATION  PROGRAM  EXPLAINED 

Available  in  Appendix  E  is  a  flow  chart  for  the  computer  simulation . • 
An  explanation  of  the  blocks  on  that  diagram  follow: 

1)  Dimension  the  computer  matrix  as  required  to  include  the  possible 
values  for  the  following  input  parameters : 
k 

XV   X2'  '"'  Xk 


(a1,  b1),  (a2,  b2),  •••,  (afc ,  bfc) 


ni'  n2'  "  *  '  ni 

2)  Generate  the  following  n.  random  repair  times  according  to  the 
gamma  distribution  with  parameters  (a.,  b.).   Utilize  the  IBM  Scientific 
Subroutine  Package  GAMMA  to  get: 


T     T     T         T 
11'   12'   13'     '   In 


1 


T     T     T     •  •  •   T 
21'   22'   23'     '   2n 


2 


T     T     T     •  •  •   T 
kl'  k2'  k3'     '  kn^. 

3)  For  each  row  above,  compute  the  mean  value  (T.)  and  the  variance 


l 
2 
(S.  )  to  get  the  pairs: 

—     2  —     2       —     2 

■p    q   .  t    q   t    q 

il'  bl  ;  i2'   2  '     '  V   k 

4)  Utilize  formula  (3-14)  to  compute  the  method  of  moments  estimate 

of  the  unknown  parameter  b ,  denoted  by  b . 

2 

5)  Compute  the  integer  closest  to  2b  for  use  with  the  X/oiTn 

distribution. 
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6)  Utilize  formula  (3-6)  to  compute  the  estimate  of  system  maintain- 
ability 0,  denoted  by  9. 

2 

7)  Utilize  formula  (3-16)  and  the  x   formula  (see  Appendix  C)  to 

compute  the  approximate  upper  confidence  limit  an.  system  maintainability. 

Denote  this  value  as  0  . 

u 

8)  Repeat  step  2  through  step  7  a  total  of  500  times. 

9)  Utilize  the  IBM  SSP  SHSORT  to  order  the  500  values  off  0  from 

u_ 

smallest  to  largest. 

10)  Pick  out  the  100th  value  of  the  above  ordering  (9  ,„0)  and  print 

ulOO 

out  this  value  as  the  20th  percentile  point  of  the  distribution.. 

11)  Compute  the  value  of  b/a  and  then  utilxze  formula  (4-3)  to  compute 
the  primary  measure  of  accuracy  for  the  simulation.   Print  out  this  value 

12)  Utilize  formula  (4-4)  to  compute  another  measure  of  simulation 
accuracy.   Print  out  this  value. 
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V.   RESULTS 

For  the  procedure  as  presented  in  Chapter  III,  the  value  of  system 
maintainability  will  follow  the  gamma  distribution.   The  parameters  of 
the  distribution  were  chosen  so  that  the  maintainability  of  each  compon- 
ent was  at  a  preselected  value  position  on  the  curve  or  would  be_  repre- 
sentative of  those  curves  currently  obtainable  in  industrial  applications, 
(See  Appendix  A.)   The  following  cases  were  simulated: 
CASE  // INPUT  PARAMETERS- 

I  k  =  15 

aJ  =  10,  b.  =1 
i        i 

X  =  .005 

n.  =  10,  20,  50,  100 
l 

II  k  -  15 

a,  =  10,  b.  =2 
i        i 

X.  =  .005 

l 

n.  =  10,  20,  50,  100 

l 

III  k  =  15 

•±  -  5,  b.  -  1 

X.  =  .005 

l 

n.  =  10,  20,  50,  100 

l 

IV  k  -  15 

a.  =30,b.  =1,  i=l,  2,  — -,10 
x         1 

a.  =  10,  b.  =  1,  i  =  11,  12,  •••,  k 

X.  =  .005 

l 

n.  =  10,  20,  50,  100 
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CASE  #  INPUT  PARAMETERS 


k  -  15 


a,  =  30,  b.  =  1,  i  -  1,  •••,  10 

a.  =  10,  b.  =  1,  i  =  11,  •••,  k 
ii 

X±   =  .01 

n.  =  10,  20,  50,  100 
l 


VI  k  =  30 


a.  =  30,  b.  =  1,  i  =  1,  •••,  10 

ii  ' 

a.  =  10,  b.  =  1,  i  =  11,  •••  ,  k 
ii  ' 

X.  =  .005 

l 

n.  =  10,  20,  50,  100 


VII  k  =  30 


a.  -  30,  b±  =  1,  i  -'  1,  •••,  10 

a.  *  10,  b.  =  1,  i  =  11,  ••-,  k 

i       i 

X.  =  .01 

1 

n.  =  10,  20,  50,  100 
l 


VIII  k  =  15 


a.  =  5,  b.  =  1 
i       l 

X.  =  .05 

l 

n.  =  10,  20,  50,  100 

l 


Each  case  is  specified  by  the  relevent  input  parameters  and  the  vary- 
ing number  of  random  repair  times  generated  (n.).   The  cases  differ 
most  significantly  in  their  values  for  the  parameters  of  the  Gamma 
distribution.   However,  the  value  of  X.  is  varied  as  is  the  value  of 
k,  the  number  of  components  in  series.   Of  course,  all  of  the  cases  were 
run  over  an  identical  range  for  n. . 
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The  results  are  compiled  in  the  following  tables.   The  column 
headings  are  self-explanatory  and  their  use  has  been  discussed  earlier 
in  this  report. 
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APPENDIX  A 


THE  GAMMA  DISTRIBUTION 
Use  of  the  gamma  function  to  describe  maintenance— time,  functions^ 
has  been  described  in  the  literature. 

The  gamma  distribution  is  a  two-parameter  function. 

b.   b.-l    x  r 
a.  i  t„  r   e 


f,  (t.;  a.,b.)  =  -   F(bJ 


The  parameters  are  denoted  by  the  letters  a  and  b.   Of  the  two,.b  is 
considered  to  be  the  more  critical  because  it  controls  the  shape  of  the 
curve;  while  a  merely  determines  the  scale  of  the  axes. 


Figure  1.   Gamma  density  function; 

a=  1,  b  =0,  1,  3,  5 

For  the  distribution  as  presented,  the  parameters  a  and  b  are- 
restricted  by  the  inequalities  a  >  0,  b  >_Q-  Also,,  both  a  and:  b:  will 
remain  as  whole  numbers. 
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The  b  parameter  is  seen  to  vary  as  a  function  of  the  repair  concept. 
When  the  repair  concept  being  described  involves  relatively  simple 
actions  (e.g.,  black  box  replacement,,  etc..)  and  very  few  long  downtime 
tasks,  the  b  parameters  will  be  a  small  number.   As  the  downtime  density 
involves  more  time-consuming  actions r   the  b  parameter  will"  take:  onr  a 
larger  value. 

Throughout  the  literature  it  has  been:  shown  that:  the  gamma:,  distribu- 
tion could  be  used  to  describe  a  family  of"  distribution  curves  varying 
from  the  negative  exponential  to  the  normal..  Thus,  by  suitable  choice 
of  the  parameters  a  and  b,  the  downtime  of  a  wide  variety  of:  reasonably 
complex  equipments  can  be  described.. 
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APPENDIX  B 

THE  EXPONENTIAL  DISTRIBUTION 
One  of  the  most  widely  used  distributions  in  fields  of  Reliability/ 
Maintainability  is  the  one-parameter  exponential  function.. 
The  function  is  defined  by: 


f(X. ,  X.)  =  X.e 
i   i     x 

=  0 


-x.x, 

x  i 


X  >  0 
elsewhere 


where  X.  is  the  failure  rate  associated  with  component  i. 


X. 


-X.t. 


-X  .X. 


r   -A.t.         —a. a. 

F(X.)  =  P(X  <  x.)  =   X.  I  e   X  X  dt.  =  1-e  X  X 
i       —  i      x  J  x 


=  0 


X.  >  0 

X- 

X.  <  0 

1_  — 


is   the  cumulative   distribution  function,  and 


E[X.]    =   1/X.. 
x  x 


Thus,  the  expected  value  equals  the  reciprocal  of  the  parameter  X 

var  (X.)  =  1/X.2 
x       x 


Figure  2.   The  exponential  distribution 
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NOTES : 

1)  The  exponential  distribution  is  a  special  case  of  both  the  gamma  and 
weibull  distributions. 

2)  It  is  the  distribution  which  is  expected  when  the  mechanisms  axe  sa 
complex  that  many  deteriorations  with  different  failure  rates  are: 
operable. 

3)  When  parts  have  an  exponential  failure  distribution,  the  equipment: 
consisting  of  these  parts  also  has  an  exponential  distribution  function. 
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APPENDIX  C 

THE  CHI-SQUARE  DISTRIBUTION 
A.  quite  useful  case  of  the  gamma  distribution  occurs  when  a  =  l/2_ 
and  b  =  n/2,  where  n  is  a  positive  integer.   Then,  a  one— parameter' 
family  of  distributions  is  obtainable  with  density  function: 


f(z)  = 


(n/2  -  1)   -z/2 
jz e 

2n/2  r  (n/2) 


z  >  0 


A  random  variable  z  having  the  above  density  is  said  to  have  the 

Z 
chi-square  distribution  with  n  degrees  of  freedom  (denoted  by  X/  \)« 

In  Figure  3  below,  the  density  function  for  n  =  1,  2,  and  n  >  2-  is 

shown. 


f(z) 


n=l 


f(z) 


fCz) 


n  >  2 


-*—  z 


Figure  3 
with  E[z]  =  n,  var  (z)  =  2n. 

Our  interest  in  the  chi-square  distribution  is  based  upon  its  many 
important  applications  in  statistical  inference. 
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The  chi-square  distribution  is  tabled  for  degrees  of  freedom  up  to 

Co.;  (a) 


30.   Thus  we  may  find  in  the  table  that  value,  denoted  by- 


satisfying  P(Z  <_  x  )  =  o,  0  <  o  <  1. 


f(z) 


X(a)(n) 


Figure  4.   The  chi-square  distribution 


For  the  case  where  the  degrees  of  freedom  exceed  30  (n  >  30)  the 
chi-square  distribution  can  accurately  be  approximated  by  the  normal 

distribution  as  indicated  in  the  following  theorem: 

2 

Theorem:      Suppose   that   the  random  variable  Y  has   distribution  x    • 


Then  for  sufficiently   large  n,    the  random  variable   /2Y  has   approximately 


the  distribution  N(/2n  -   1,    1).      Therefore, 


26 


P(Y  1  t)   =  P(/2Y  £  *^2t) 

=  P(/2Y  -   /2n  -  T  <_  Tit  -   /2n  -  I) 
=  $(/2t"  -   /2n  -  1) 
where   $  values   are   obtained   from  the  normal   tables. 
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APPENDIX  D 

THE  METHOD  OF  MOMENTS 

The  oldest  general  method  for  generating  estimates  of  unknown: 
parameters,  given  a  sample,  is  known  as  the  method  of  moments.   It,:  is 
a  comparatively  simple  method  and  it  generates  "reasonable"  estimators. 

In  general,  the  method  proceeds  as  follows: 
Given  a  random  variable  X,  with  distribution  function  F  wher&  this 

X 

distribution  is  indexed  by  the  unknown  parameter  X.   Assume  the  first 

moment  of  X  (its  mean  value)  is  dependent  in  some  simple  way  upon  X, 

like  uv  =  g(X).   Then,  given  a  sample  of  n  values  of  X,  define  the  first 
x 

sample  moment  as  X.   The  method  of  moments  now  specifies  that  X  be  equated 
to  g(X)  and  finally  to  solve  for  X.   This  resulting  value  for  X  is  the 
method  of  moments  estimator  of  X. 

For  the  two  dimensional  case,  such  as  X  being  distributed  according 
to  the  gamma  distribution  with  unknown  parameters  a  and  8,  there  is 
required  a  multi-dimensional  parameter  space  Q  =   E  which  in  the  gamma 
case  has  k  =  2.   Then,  this  space  can  be  defined  such  that 
ft  =  {(a,  6);  a  >  0,  8  >  0}  «3  e  and  further,  it  can  be  shown  that 

2     2    2  2 

u  =  a8  and  y   =  a0   +  a  8 
X  X 

2 
where  y    is  the  second  moment  with  respect  to  X.   Then,  as  in  the  one 

A 

dimensional  case,  a  random  sample  is  required  which  can  be  analyzed  to 
give: 

X  =  ct8 

2—2     2    2  2 

S^  +  X  =  aB  +  a  0 


28 


which  allows ,  as  shown  in  Chapter  III 

-2  S2 

sx  x 

as  acceptable  estimators  of  the  parameters  a  and 
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APPENDIX  E 


SUBROUTINE 
GAMMA 


® 


READ  IN 
PARAMETERS 


GENERATE 

GAMMA 

VARIATES 


COMPUTE 
T  and  S' 


COMPUTE 
b  and  2b 


COMPUTE 
9 


UTILIZE  x 
TO  COMPUTE 


SUBROUTINE  V 
SHSORT     / 


ORDER  500 
VALUES  OF 

u(i) 


PRINT  20th 
PERCENTILE 


COMPUTE 
b/a 


COMPUTE 

[#e  ...  >  b/a] 
u(i)  ~ 

500 


PRINT  OUT 
STATISTICS 


(      STOP  J 


© 
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C      THIS  PROGRAM  COMPUTES  THE  80  PERCENT  UPPER  CONFIDENCE  LIMIT 
C      DIMENSION  THE  MATRIX 

REAL*4  LAM6DM  100) ,N( 100)  ,LAM 

DIMENSION   T(  100, 100)    ,TB ( 100 ) , S2( ICO >  ,TU(5C0)  ,ALPHA(  100)  ,  BETA (10 
$0),KEY(500) 
C      READ  INTO  PROGRAM  PARAMETERS  OP  INTEREST 
12  35  READ(5t  lOOf  END=4321  )  K,  (LAMBDA(  I )  ,I=1,K) ,  (ALPHA(J)  r  BETA  (  J  )  ,J  =  1 ,  K  ) 
READ(5,105)  MI 

PROGRAM  WILL  LOOP  500  TIMES,  GENERATING  INDEPENDENT  ESTIMATES  OF 
THE  80  PERCENT  UPPER  CONFIDENCE  LIMIT 

DO  1234  JJ=1,500 

GENERATE  GAMMA  TIMES  TO  REPAIR 

DO  1  1=1, K 

N( I)=NI 

NN=NU  ) 

DO  1  J  =  1,NN 

CALL  GAMMA (BET  A (  I  )  ,  ALPHA  (  I)  ,X  ) 

T(I ,J)  =  X 

1  CONTINUE 

COMPUTE  MEAN  TIME  TO  REPAIR  FHR  EACH  K  GAMMAS  GENERATED  AND  THE 
VARIANCE  FROM  THE  METHOD  OF  MOMENTS  FIT 

DO  2  1=1, K 

SUM=0.0 

SUMSO=0.0 

NN=N(I ) 

DO  3  J=1,NN 

SUM=SUM  +  T( I, J) 

3  SUMSQ  =  SUMSQ  +  T ( I , J ) *T (  I , J  ) 
TB(  I  )  =  SUM/M(  I ) 
S2(I)  =  (SUMSQ-N( I )*TB(  I)*TB(  I)  )/(N( I )-1.0) 

2  CONTINUE 
LAM=0.0 
DO  4  1=1, K 
LAM=  LAM+LAMBDA( I) 

4  CONTINUE 

SUM=0.0 

SUMSQ=0.0 

DO  5  1=1  ,K 

SUM=  SUM  +  (LAMBDA( I) /LAM)*TB( I ) 

SUMSO  =  SUMSQ  +  (LAMBDA! I )/LAM)*(LAMBDA(I)/LAM)*(S2(I)*S2(I ) ) 
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5  CONTINUE 
C 

C      COMPUTE  VALUE  FCP  ESTIMATE  OF  PARAMETER  BETA,  CALLED  BHAT 
C 

BHAT=(SUM*SUM) /SUMSQ 
C 

C      COMPUTE  INTEGER  CLOSEST  TO  TWO  BHAT 
C 

IH=2.0*BHAT+0.5 

BHAT2  =  IH 
C 

C      COMPUTE  ESTIMATE  OF  MEAN  TIME  TO  REPAIR  FOR  SVSTEM,CALL  IT  THETA 
C 

THETA  =  0.0 

DO  6  1=1 ,K 

THETA  =  THETA  +  ( L AMBDA( I )/L AM >*T8f I ) 

6  CONTINUE 
C 

C      UTILIZE  CHISQ  FORMULATION  TO  GET  STATISTIC  FOR  CONFIDENCE  LIMITS 
C 

CHISA=(2.0/(9.0*BHAT2)  ) 

CH I SQ=  B  HAT2* ( 1 . O-CH I S A-C . 341 78*SQRT ( CH I S A ) ) **3 

TU( JJ)=( THETA* BHAT 2) /CHISQ 
C 
C 

1234  CONTINUE 
C 

C      NOW  THAT  THE  500  ESTIMATED  VALUES  HAVE  BFEN  GENERATED,  SORT  THEM 
C      INTO  ORDER 
C 

KK=500 

DO  9<*  I=1,KK 
99  KEY(I)=I 

CALL  SHSORT(TU,KEY,KK) 
C 

C      PICK  OUT  AND  PRINT  THE  2CTH  PERCENTILE  ACCORDING  TO  THE  FOLLOW- 
C      ING  FORMAT 
C 

WRITE(6,?22)  M I , K, LAMBDA ( 1 ) 

WRITE(6,200)  TU(IOO) 


C 

C 


222  FORMAT(  ■  l1  ,  «NI  =  M5,  «  K=*,I5,  •  LAM6DA=  *  ,  F10.  5,  ///  ) 

100  F0RMATU5/  15F5.3  /  15F5.3  /( 10<P4.0, F4. 0 ) ) ) 

105  FORMAT ( 15) 

200  FORMAT  (IX,'  THE  20TH  PERCENTILE  IS  »,F15.5) 

WRITE ( 6,400) <TU( I ) , 1  =  1  ,500) 
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400  FORMATS  •  ,  10(  F8.4  ,  2X  )  ) 
C 

C      COMPUTE  THE  VALUE  OF  THE  THEORETICAL  MEAN,  CALL  IT  BOA 
C 

DENOM=0.0 

TOP=0.0 

DO  50  1=1, K 

TOP  =  TOP   +(LAMBDA(I)/LAM)*(BETA(I )/ALPHA(I) ) 

DENOM=DENOM+(L A^BDA( I) /LAM )**2  *< BETA ( I ) / ALPH A{ I ) **2) 
50  CONTINUE 

A  =  TOP /DEMON 

B  =  T0P*-*2/DEN0M 

BOA  =  B/A 
C      COMPUTE  THE  ACTUAL  PERCENTAGE 

C=0.0 

DO  60  1=1,500 

IF  (TU(  I  ).GE.BOA)  C=C+1.0 
60  CONTINUE 

STAT  =  C/500.0 
C      COMPUTE  THE  ABSOLUTE  VALUE  OF  THE  RELATIVE  DIFFERENCE 

E=   ABS(TU(100)-80A)/B0A 
C      PRINT  OUT  THE  ABOVE  STATISTICS  ,  ACCORDING  TO  THE  FOLLOWING  FORMAT 

WRITE(6,?20)  BOA, STAT, 3, A 
220  FORMATdX,'  B  OVER  A  IS  ',F15.5,5X,'  STAT  IS  «  ,F1 5. 5,5X  ,••  B  IS  «, 
$F15.5,5X, •  A  IS  ■ ,F15.5) 

WRITE(6,230)  E 
230  F0RMAT(1X,«  ABS  OF  VALUE  IS  «,F15.5) 

IN  =  20 

RANGE(l)  =   0.1822  +  0.00159 

DO  91  1=2,20 
91  RANGE(I)  =  RANGE(I-l)  +  0.00159 

K=l 

DO  98  1=1,20 

FREQ( I )  =  0.0 

DO  97  J=K,500 

IFCRANGEU  ).LT.TU(  J)  )  GO  TO  96 

FREQf I )=FREO( I J+1.0 

GO  TO  97 

96  K=J 

GO  TO  9  8 

97  CONTINUE 

98  CONTINUE 
GO  TO  1235 

4321  STOP 
END 
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SUBROUTINE    GAMMA(B,AtX) 

K=B 

TR=1.0 

DO    5    I=1,K 

R=RN(0) 
5    TR=TR*P 

X=-ALQG(TR)/A 

RETURN 

END 

******************************************************* ***************^ 
afc 

*  PURPQSF 

*  RN    COMPUTES    UNIFORMLY    DISTRIBUTED    RANDOM    VARIATES    OVER    THP    RANGE 

*  (0,1).       THE     PERIOD    FOR    THIS    GENERATOR     IS     2**32    WHICH    IS    THE    LENGTH 

*  OF    THE    FULL     SYSTEM/360    WORD    GIVING    TWICE    THE    PERIOD    OF    RANDU. 

*  USAGE 

*  IN  X    =    RN(O)        WHERE    X    WILL    BE    ASSIGNED    THE    NEXT    RANDOM    DEVIATF. 
* 

*  DESCRIPTION    0^    PAPAMETFRS 

*  NO    FORMAL    PARAMETERS 
* 

*  REMARKS 

*  THIS    ASSEMBLY    LANGUAGE    RnUTINE    GIVES    TWICE    THE    PERIOD    OF     EXISTING 

*  GENERATORS    AND     TS    EXTREMELY    FAST.        11.88    MICROSECONDS    WITH    4.23 

*  MICROSECONDS    FOR     TH^    CALLING    SEQUENCE.        IT     IS    TAKEN    DIRECTLY    FROM 

*  ,     COMM.     ACM     12,12     (DEC.     1969)     695. 
* 

*  METHOD 

*  LEHMER'S    MULTIPLICATIVE    CONGRUENTIAL    SCHEME     IS    USED. 

*  U(K+1)     =     L*U(K)      (MOD    P),        X(K+1)     =    U(K+1)/P       WHERE    P    =    2**32    AND 

*  L    =    32781. 
* 
RN 

SET    BASE     ADDRESS 

LOAD    U(K) 

SEED  =  SEEO*M'JLT 

STORE    U(K+1) 

FLOAT    AND 

NORMALIZE 

RETURN 
MULT  DC  F' 32781 «        2**32 

ZERO 
FLOT 
SEED 
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CSECT 

USING 

*,15 

L 

1,SEED 

M 

0,MULT 

ST 

I, SEED 

LD 

0,FLOT 

AD 

0,ZE«0 

BR 

14 

DC 

F'327811 

DC 

D'O' 

DC 

X'46000000 

DC 

F«  1» 

END 

RN 

RN 
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